Influence of the Potential Energy Landscape on the Equihbration and Speciflc Heat of 

Glass Forming Liquids 
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We show that a glass transition, signaled by a peak in the specific heat vs. temperature, can 
occur because a glassy system that shows no signs of aging progresses so slowly through the energy 
landscape that the time needed to obtain an accurate estimate of the thermodynamic averages 
exceeds the observation time. We find that below the glass transition temperature of a three 
dimensional binary mixture of soft spheres, the specific heat increases with measurement time spans 
orders of magnitude longer than previously recognized equilibration times. 



As a glass forming liquid is cooled or compressed, a 
glass transition occurs when the system falls out of equi- 
librium, i.e., when the time scale for reaching equilibrium 
exceeds the observation time 0| . Like many complex sys- 
tems, such as proteins and neural networks, the dynamics 
of such a system is strongly influenced by the potential 
energy landscape where each point corresponds to a par- 
ticular configuration and energy of the system ^ . The 
energy landscape can be used to describe the three ways 
in which a system can fall out of equilibrium. First a sys- 
tem can become trapped in a mctastable minimum where 
it stays for the duration of the observation time. Second 
a system can be in an energetically unlikely part of phase 
space and proceed slowly to a region of the energy land- 
scape where its configurations obey a Boltzmann distri- 
bution. As such a nonequilibrium system evolves toward 
more probable regions of phase space, it exhibits aging 
which means its properties systematically change with 
time and do not obey stationarity Q. The aging time, 
after which aging stops, is equal to the a relaxation time 
which is the charactistic time for the system to forget its 
initial configuration. 

The third way is not widely appreciated and is the sub- 
ject of this paper. Namely, even after a glassy system no 
longer ages and has reached basins with appropriate en- 
ergies, the system proceeds so slowly through the energy 
landscape that it takes a long time to accumulate the 
large number of statistically independent measurements 
needed to accurately determine a thermodynamic aver- 
age. We find from molecular dynamics simulations that 
a glass forming liquid can undergo a glass transition, as 
signaled by a peak in the specific heat Cy versus temper- 
ature, that is due to insufficient averaging by a system 
that shows no signs of aging. The distribution of energies 
that a system samples in a basin of the energy landscape 
is a subset of the full distribution of energies available to 
the system. Since this subset has a smaller variance than 
the full distribution, the resulting specific heat, which is 
proportional to the variance of the energy, will be smaller 
when calculated from short time spans than from long 
time spans. These smaller values account for the values 
below the peak in Cy on the low temperature side. Going 



to longer time spans eliminates the peak, though at tem- 
peratures below the peak temperature Tp, these longer 
time spans can be orders of magnitude longer than pre- 
viously recognized equilibration times such as the a re- 
laxation time, the energy correlation time, and the aging 
time. 

We have performed a molecular dynamics simulation 
on a three dimensional glass forming liquid |, I con- 
sisting of a 50:50 mixture of two types of soft spheres, 
labelled A and B, which differ only in their sizes. The 
interaction between two particles a distance r apart is 
given by Vapir) = e[{cj ap / rf"^ + ^apir)] where the in- 
teraction length Gap = (ca + o'/3)/2, and the ratio of the 
ratio of the diameters (JbIcta = 1-4 {a, (3 = A, B). The 
cutoff function X^fiir) = r/aa/s - A with A = 13/12i2/i3, 
The interaction is cutoff at the minimum of the poten- 
tial Vaf}{r). Energy and length are measured in units of 
e and a a, respectively. Temperature is in units of e/fc^, 
and time is in units of a a y^m/e where m, the mass of the 
particles, is set to unity. The equations of motion were 
integrated using the leapfrog method |Q with a time step 
of 0.005. During each run the density Po = N/L^ was 
kept constant. N = Na + Nb is the total number of 
particles. The system occupies a cube with volume L'^ 
and periodic boundary conditions. According to the ideal 
mode coupling theory, the relaxation time diverges at a 
temperature Tc (|]. For our system Tc = 0.303 §. 

We cool the system from a high temperature (T=:1.5) 
by lowering the temperature in steps of AT ~ 0.05. At 
each temperature we equilibrate for 10** time steps and 
then measure the quantities of interest for 10^ additional 
steps. The results are then averaged over different runs. 
At the glass transition the system falls out of equilibrium 
and becomes trapped in a basin of the energy landscape. 
In order to avoid this, we have used parallel tempering 
together with molecular dynamics. We implement paral- 
lel tempering (FT) |9[ ^ by running molecular dynam- 
ics simulations in parallel at chosen temperatures using 
a temperature constraint algorithm to keep the tem- 
perature of each simulation constant. At 100 time step 
intervals we attempt to switch the configurations of two 
neighboring temperatures using a Metropolis test which 
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ensures that the energies of the configurations sampled 
at any given temperature have a Boltzmann distribution. 
Let Pi and P2 be two neigboring inverse temperatures, 
and let Ui and U2 be the corresponding potential en- 
ergies of the configurations at these temperatures. If 
A = (/?! — fi2)[U2 ~ Ui), then the switch is accepted 
with probability unity if A < and with probability 
exp(— A) if A > 0. The acceptance ratio is between 30% 
and 75%. Near Tp, the acceptance ratio was above 60%. 
After a swap is accepted, the velocities of the particles in 
each configuration are rescaled to suit their new temper- 
ature. Switching configurations allows a given simulation 
to do a random walk in temperature space in which it vis- 
its both low temperatures and high temperatures. This 
helps to prevent it from becoming trapped in a valley of 
the energy landscape at low temperatures. Typically we 
equilibrate for 2 million time steps and then do measure- 
ments for an additional 4 to 10 million molecular dynam- 
ics (md) steps. We average over different runs which have 
different initial positions and velocities of the particles at 
each temperature. 

We calculate the specific heat Cy per particle at con- 
stant volume in two ways. The first uses the fluctu- 
ations in the potential energy U per particle: Cy = 
(3A:b/2) + NkB0^ {{U^} - {Uf) where the first term is 
the kinetic energy. The second way uses Cy = d < U > 
/dT « (< f/(r.+i) > - < U{T,) >)/(T,+i -T,). The re- 
sults are shown in Figure |l]. Notice that there is a sharp 
asymmetric peak centered at Tp = 0.305 ±0.003. The low 
temperature side of the peak drops steeply. The curves 
for 216 and 512 particles coincide, indicating that there 
is no size dependence. The discrete points are calculated 
from the energy fluctuations. The solid line is calculated 
from the derivative of the energy. The fact that the two 
coincide indicates that the system was equilibrated in all 
the basins of the energy landscape that were visited. For 
comparison we also show the result of cooling through 
the transition (o). The peak in C'v found by cooling 
coincides with the peak found in parallel tempering. 

We do not believe that the parallel tempering peak in 
the specific heat is due to the system becoming trapped 
in a metastable minimum in the energy landscape for 
the following reason. The configurations corresponding 
to the bottom of the basins of the energy landscape are 
called inherent structures ||ll[. We sampled the config- 
urations that were visited during the parallel tempering 
runs and found the corresponding inherent structure en- 
ergy e/s per particle by minimizing the potential energy 
locally using the method of conjugate gradients ||l^. In 
the inset of Figure |l| we plot the average inherent struc- 
ture energies versus the temperature of the configuration 
that was originally saved. As the temperature decreases 
below 0.5, (e/s) decreases rather steeply 0. The tem- 
perature of the inflection point of this decrease coincides 
with the temperature Tp of the peak in the specific heat. 
We also show (e/s) for a system of 512 particles that 
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FIG. 1: Specific heat versus temperature. Shown is the spe- 
cific heat calculated from fiuctuations for systems with 216 
particles (Q, averaged over 6 runs) and with 512 particles for 
measurements covering 4 x 10® time steps (v, averaged over 9 
runs) and 10^ time steps (□, averaged over 3 runs) done with 
parallel tempering. The solid line is the specific heat calcu- 
lated from the derivative of the energy from the 4 million time 
step parallel tempering runs with 512 particles. Also shown is 
the specific heat calculated from the energy fluctuations for a 
system with 512 particles (o, averaged over 6 runs) that was 
cooled conventionally with 10^ time steps per temperature. 
The ■< correspond to concatenating between 13 and 40 single 
temperature runs; each run had 10* time steps with 512 parti- 
cles that were initiated from parallel tempering configurations 
and equilibrated for up to 5 x 10® time steps. Inset: Average 
inherent structure energy per particle versus temperature for 
512 particles obtained from parallel tempering, cooling, and 
single temperature runs. The symbols denote the same cases 
as in the main figure. The arrow points to the inflection point 
which coincides with Tp. 



was cooled from T = 1.5. At low temperatures (e/5) 
is rather independent of temperature for the cooled sys- 
tem, indicating that the system is trapped in an energy 
basin. Such a flattening off at low temperatures is not ob- 
served when parallel tempering is used, indicating that 
the system is able to visit deeper energy basins as the 
temperature decreases. 

However, this does not mean that the peak in the spe- 
cific heat is an equilibrium feature. We will now show 
that the peak is the result of not sampling enough of 
phase space below Tp. None of the runs we have done at 
T < Tp have reached the time which is necessary in order 
to achieve an accurate thermodynamic average. To show 
this, we took a configuration of 512 particles generated 
by parallel tempering at T = 0.289855, equilibrated for 
5 X 10^ time steps, and then ran for an additional 10* time 
steps during which we recorded the energy at every time 
step. Then we did block averaging in which we divided 
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FIG. 2: Block averaged specific heat versus time span Att 
for 512 particles at T = 0.289855, 0.308642, and 0.3424658. 
Solid lines are the result of block averaging each run and then 
averaging over the number of runs shown. Open symbols are 
the result of stringing runs together and then block averaging. 
Inset:Block averaged specific heat versus time for 512 particles 
at T = 0.289855. The lower curve corresponds to Ah = 10** 
md steps, and the upper to Att = 4 x 10® md steps. The time 
is the time (in md steps) in the middle of each block. The 
data at T = 0.289855 is averaged over 23 runs. Parameters 
for both figures are the same as in Figure 0. 



our 100 million time steps into equal segments, each of 
length Atb, and calculated the specific heat from energy 
fluctuations for each segment |l3). The block averaged 
specific heat versus time for 2 different time spans Atb is 
shown in the inset of Figure |^. Notice that for any given 
time span, there is no sign of systematic aging. However, 
the specific heat averaged over time increases with A^f,. 

The specific heat, averaged over time spans of a given 
size and over different runs, versus time span size Atb at 
several temperatures is shown in Figure || by the solid 
lines. To obtain time spans that are longer than any 
given run, we concatenated the energies from the runs 
done at a given temperature to make one huge run, and 
then did block averaging on the huge run. The results 
are shown as open symbols in Figure |^. There is good 
agreement between the solid lines and symbols. Notice 
that the specific heat initially increases with time span 
but then levels off when the time span is long enough to 
exceed the time needed to achieve the thermodynamic 
value of the specific heat. This time increases with de- 
creasing temperature. Thus at T = 0.289855 < Tp, the 
specific heat continues to increase with time span up to 
Atb = 10^ time steps which is 100 times longer than the 
a relaxation time t. Equilibrated values of Cy near Tp 
are plotted in Figure ^ and lie above the peak in Cy 
found with parallel tempering. Thus the specific heat 



peak found with parallel tempering is the result of not 
sampling enough of phase space at T < Tp to obtain the 
true thermodynamic value of the specific heat Cy^'^. The 
exploration of the energy distribution below Tp is slow 
even with parallel tempering because the probability of 
sampling large increases in U are exponentially small. 

Note that Cy"° is proportional to the variance afj of 
the distribution P(C/), i.e., C^J"<= = Nalj/ikeT^)- If cr^ 
is finite and if n sample values of U are statistically inde- 
pendent and identically distributed, then basic statistics 
dictates that the measured Cy, which is proportional to 
the sample variance of U, has an expectation value of 
{Cv) = C^™'=(1 - 1/n) |l|. Since the potential energies 
in a time series can be correlated, the number n of sta- 
tistically independent energies is given by n — Atb/ru 
where tu is the energy correlation time. Fitting the 
data that is within 5 to 10% of Cy"^ in Figure | to 
{Cv) = C^™<=(1 - l/n) yields tc/ « 3 x 10^ 1 x 10^, 
and 5 X 103 time steps at T = 0.2898551, 0.308642, and 
0.342466, respectively. These values are comparable to 
the a relaxation times r of 1 x 10^, 6 x lO'^, and 2 x lO'^ 
time steps, respectively (see below). They imply that the 
energies are correlated and that n is substantially smaller 
than the total number of energies. 

Note that at shorter time spans and lower tempera- 
tures Figure |2| shows that Cy ~ In(Atb). This occurs be- 
cause the system does not uniformly sample P{U) during 
these shorter time spans. As the system travels through 
the energy landscape, it samples the energies of each 
basin that it visits. As shown in Figure ^, we find that the 
distribution of energies sampled from the basins visited 
during a shorter time span has a smaller variance than 
the total distribution P{U). Furthermore, the centers of 
the smaller distributions do not necessarily coincide with 
the center of the total distribution. Rather the smaller 
distributions are centered at the inherent structure en- 
ergy plus the energy of vibrations around ejs |jl^. In 
support of this, we show in the inset of Figure ^ that the 
block averaged energy versus time approximately coin- 
cides with {eis{t) + 3A:bT/2) vs. time. As more basins 
are visited, the sample average (U) moves towards the 
average of the full distribution and the variance grows. 
This corresponds to Cy increasing with time span. 

We now cite evidence that systems at temperatures 
just below Tp have equilibrated in the sense of showing 
no signs of aging. First the inset of Fig. ^ shows the lack 
of aging in the specific heat versus time for a given value 
of Atb- Second is the absence of aging in the a relaxation 
time T. We have calculated r using the full intermediate 
scattering function FBsik^t) — {1 /Nb) {pj:{t) p _^{0)) for 
the B particles where Nb is the number of B particles and 
Pfc(*) = S£i exp(-zA- ■ ri{t)). We choose k = /cmax, the 
wavevector of the maximum in the partial static struc- 
ture factor SBB{k) for the B particles, because FBB{k, t) 
relaxes slowest at fcmax |15|. The relaxation time r is de- 
fined by F{k,T)/ F{k,t = 0) = 1/e. We have averaged 



4 



5e+06 



4e+06 



q 
■(5 

^ 3e+06 

c 

o 

O 

O 2e+06 
E 

^ 1e+06 









2e+07 4e+07 6e+07 8e+07 1e+08 

Time [md steps] 



1.5 1.6 1.7 1.8 1.9 

Average Potential Energy Per Particle 



FIG. 3: Number of configurations versus tlie average energy 
per particle acquired during one run of 10** md steps of a sys- 
tem of 512 particles at T=0. 289855. The large curve centered 
at 1.62 is the histogram for the entire run. The small curve 
on the left centered at 1.60 represents the counts acquired 
between 4 x 10^ and 5 x lO'^ time steps; the small curve on 
the right centered at 1.63 represents the counts acquired be- 
tween 7 X lO'^ and 8 x lO'^ time steps. Inset; Data from the 
same run. Circles represent the potential energy per particle 
averaged over lO'' time steps versus time. The solid line is 
eis + 3fcsr/2 versus time. 



F{k,t)/F{k,t = 0) over 40 runs after waiting times tw 
of 5 X 10^, 10^, and 1.5 x 10® time steps for a system 
with 512 particles at T = 0.289855 < §. We find 
that the a relaxation time r = (1.0 ± 0.1) x 10® md time 
steps. This value of r shows no signs of aging Q in the 
sense that there is no systematic variation with waiting 
time. The lack of aging is to be expected since the aging 
time is equal to r which is much less than tw- We have 
confirmed that the aging time is the same as the equilib- 
rium value of T by starting from 11 different equilibrium 
configurations at T = 1.5, quenching to T = 0.289855, 
and measuring r after waiting times of tw = 0, 5 x 10'*, 
5 X 10^ 10^ 5 X 10^ and 2.5 x lO"^ md steps. For tw < 10® 
md steps, r increases with tw which indicates aging [Q. 
However, for tw > 10® md steps, there is no aging and r 
equals its equilibrium value of 10® md steps. 

Third we have looked for signs of aging in the in- 
herent structure energy versus time in our single tem- 
perature molecular dynamics runs of 10® time steps at 
T = 0.289855 < Tp. During the run, configurations were 
recorded every so often. Averaging over 40 runs, we find 
no evidence that the inherent structure energy decreases 
systematically with time, though the noise in the data 
prevents us from seeing changes smaller than 1%. 

We have also examined the root mean square displace- 
ment (Ar2(t))i/2 ^here {Ar^it)) = (l/^)(Eti(r^(i) " 



examined after 10® time steps at T = 0.289855. This is 
comparable to the box size L = 9.48(7^ for a system of 512 
particles. So the system does not appear to be getting 
stuck in a metastable minimum of the energy landscape. 

Our work is a cautionary tale for those who perform 
numerical simulations on slowly relaxing systems. It in- 
dicates that to obtain accurate thermodynamic averages, 
one must not only check that the system shows no signs 
of aging, but one must also check that the quantity to be 
measured has sampled enough of phase space to obtain 
a large number of statistically independent values. This 
sampling time can be orders of magnitude longer than 
previously recognized time scales such as the aging time 
and the a relaxation time. 
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